This lecture will walk you through some more advanced spatial data manipulation and visualization in R. We acknowledge that the materials provided here do not by any means intend to be comprehensive nor attempt to cover all concepts that would otherwise require a semester-long course. The sections and code chunks found here are derived from a variety of internet resources as well as past R workshops. Because this repository is open-source, please feel free to use it as you deem appropriate and I would greatly appreciate any feedback to improve and build upon it. If you are familiar with version control concepts, such as Git and Github, you are more than welcome to ‘fork’ this repository on your own account and send a ‘pull request’ to me if you wish to add/edit any of the scripts.
Here, we’ll use the setwd function to tell R where the module 1 folder is located on our machine. This will depend on the location of your downloaded course folder, Intro-Spatial-R.
# modify the filepath based on your user specific folder location
# Single or double quotes work for text strings. Must use forward slashes.
#Some common examples:
#setwd('C:/Users/username/Documents/GIS610_SpatialR/Spatial Data Visualization in R') # for Windows users
#setwd("~/Documents/GIS610_SpatialR/Spatial Data Visualization in R") # ~ for Mac users
Make sure you have these packages installed (use install.packages() function in R) or some scripts will not run. In order to use some of the useful functionalities that come with each package, we need to call them into our current R session using the library() or require() command.
library(rgdal) #for reading/writing geo files
library(rgeos) #for simplification
library(plyr)
library(leafletR)
library(RColorBrewer)
library(Quandl)
library(reshape2)
library(RgoogleMaps)
library(plotGoogleMaps)
library(loa)
library(latticeExtra)
library(ggmap)
library(maptools)
library(rjson)
library(googleVis)
This lecture will also touch on topics such as interactive web maps from R. This topic is fairly recent in the R community and package development is ongoing. In fact, you will not be able to find and install the following packages from the common CRAN repository. In order to download and install R packages that are under development on the author’s Github repository, we need to install an R library called devtools, which allows to use the command install_github() for the desired packages.
library(devtools)
## Warning: package 'devtools' was built under R version 3.2.3
##The R packages rCharts and rMaps are not available on CRAN yet.
##Run the following code to install from Github repo prior to loading the libraries:
##install_github('ramnathv/rCharts')
##install_github('ramnathv/rMaps')
library(rCharts)
library(rMaps)
In the previous lecture, we saw examples on how to produce basic maps within R as well as adding some more advanced descriptive information to them. In this section, we keep showing alternative ways in which R can be used as a GIS, i.e. for basic-to-intermediate mapping tasks. Let us start with some examples on how Google Maps basemaps can be pulled into R and basic maps be overlaid onto these “canvases”. Two very useful R packages are RgoogleMaps and loa.
RgoogleMaps and loa packagesThe three core functions within RgoogleMap are:
PBSmapping to an existing map in form.#load dataset of interest with crime data in the SF area for 2013
load("./data/incidents2013.rda")
#list what dataset(s) has been loaded
ls()
## [1] "incidents"
#let's have a quick exploration of the crime dataset
dim(incidents)
## [1] 126991 17
head(incidents)
## IncidntNum Category Descript DayOfWeek
## 1 121047551 RECOVERED VEHICLE VEHICLE, RECOVERED, AUTO Saturday
## 2 13670901 LARCENY/THEFT PETTY THEFT OF PROPERTY Monday
## 3 60847080 WARRANTS WARRANT ARREST Sunday
## 4 71276391 WARRANTS ENROUTE TO OUTSIDE JURISDICTION Wednesday
## 5 71276391 WARRANTS WARRANT ARREST Wednesday
## 6 90818299 WEAPON LAWS POSS OF LOADED FIREARM Wednesday
## Date Time PdDistrict Resolution Location
## 1 2013-01-05 12:18 BAYVIEW ARREST, BOOKED 900 Block of CONNECTICUT ST
## 2 2013-04-15 16:15 INGLESIDE NONE SILVER AV / MISSION ST
## 3 2013-03-24 17:33 SOUTHERN ARREST, BOOKED 800 Block of BRYANT ST
## 4 2013-02-13 22:41 NORTHERN ARREST, BOOKED CHESTNUT ST / FILLMORE ST
## 5 2013-02-13 22:41 NORTHERN ARREST, BOOKED CHESTNUT ST / FILLMORE ST
## 6 2013-02-20 16:00 SOUTHERN ARREST, BOOKED 800 Block of BRYANT ST
## X Y violent HrOfDay TimeOfDay HourOfWeek PhaseOfWeek
## 1 -122.3978 37.75397 FALSE 12 12.30000 36.3000 Saturday
## 2 -122.4313 37.72873 FALSE 16 16.25000 88.2500 Mon-Thu
## 3 -122.4037 37.77518 FALSE 17 17.55000 65.5500 Sunday
## 4 -122.4363 37.80081 FALSE 22 22.68333 142.6833 Mon-Thu
## 5 -122.4363 37.80081 FALSE 22 22.68333 142.6833 Mon-Thu
## 6 -122.4037 37.77518 FALSE 16 16.00000 136.0000 Mon-Thu
## CensusBlock
## 1 06075061400
## 2 06075025500
## 3 06075018000
## 4 06075012800
## 5 06075012800
## 6 06075018000
#let's fetch the map tile at the desired resolution and centered on a lat/lon point
SFzoom13 <- GetMap(center = c(lat = 37.77173, lon = -122.4306),
destfile = "./data/figures/SanFrancisco.z13.png",
GRAYSCALE = FALSE, zoom = 13, SCALE = 2, maptype = "terrain")
#let's randomly sample a number of crimes for visualization purposes
Ncrimes <- 10000
ranPts <- sample(nrow(incidents), Ncrimes)
crimes <- incidents[ranPts, ]
#We are mostly interested in looking at "violent" crimes.
#let's start by converting the variable from factor to logical
crimes$violent <- as.logical(crimes$violent)
#We want to pick two different colors for violent and non violent crimes.
#in order to overlay these colors with some transparency, we need to define the transparency "alpha" level
#inside the rgb() function. In this case we pick the colors 'red' and 'cyan' defined by
#their corresponding R-G-B values
colViolent <- rgb(1,0,0,0.3)
colNonViolent <- rgb(0,0,1,0.3)
cols <- c(colViolent, colNonViolent)
#now, we are ready to plot on our static Google Maps tile
#png("./data/figures/SanFrancisco.z13.png", width=1280, height=1280) #use if you want to save higher quality image
tmp <- PlotOnStaticMap(SFzoom13, lat = crimes$Y, lon = crimes$X, cex=0.9,
pch=20, col=cols[as.numeric(crimes$violent)+1])
legend("topright", col = c(colNonViolent,colViolent), pch=20,
legend = c("N","Y"), title="violent", cex = 1)
#dev.off() #if you use png() remember to close the graphical device
The loa package contains various plots and functions that make use of the lattice/trellis plotting framework. The plots (which include loaPlot, GoogleMap and trianglePlot) use panelPal, a function that extends lattice and hexbin package methods to automate plot subscripting and panel-to-panel and panel-to key synchronization/management. See ?loa for further details. Here, we will use GoogleMap(), which is basically a wrapper for RgoogleMaps function GetMap() and loa function loaPlot()
#let's save the time of the day in a separate variable and look at it
temp <- incidents$TimeOfDay
head(temp)
## [1] 12.30000 16.25000 17.55000 22.68333 22.68333 16.00000
#we not want to generate a character vector with labels indicating the approx. time of the day
#We can start by replicating the label "About Midnight" for a number of times equal to the length
#of the temp vector
tod <- rep("About Midnight", length(temp))
#now that we have prepared this new character vector variable, we need to substitute
#the correct labels to the appropriate time windows in a day. We can do it in the following way
tod <- ifelse(as.numeric(temp) > 2 & as.numeric(temp) < 9, "About Dawn", tod)
tod <- ifelse(as.numeric(temp) > 8 & as.numeric(temp) < 14, "About Midday", tod)
tod <- ifelse(as.numeric(temp) > 13 & as.numeric(temp) < 20, "About Dusk", tod)
tod <- factor(tod, levels=c("About Dawn", "About Midday", "About Dusk", "About Midnight"))
#we want to also have a character vector with a label indicating whether a crime is violent or not
v <- ifelse(incidents$violent==TRUE, "Violent", "Non-violent")
#now we are redy to use the GoogleMap() function to plot a series of temporal plots
#of crimes by time of the day and type of crime (violent/non-violent)
GoogleMap(~Y*X|tod+v, #conditioning
data=incidents, col="darkred",
cex=0.1, alpha=0.1, pch = '.',
scales=list(draw=FALSE), xlab="", ylab="") #to hide axes
#png("SFtest.png", 1690, 717, bg="transparent", res=118, type="cairo-png")
useOuterStrips(trellis.last.object())
#dev.off()
ggmap packageThe ggmap package enables visualization by combining the spatial information of static maps from Google Maps, OpenStreetMap, Stamen Maps or CloudMade Maps with the layered grammar of graphics implementation of ggplot2. More information on this R package and tutorials from which the following section is extracted can be found HERE.
#let's start simple by using the ggmap() function to simply pull down a map from the web
#centered on a city of interest
raleigh <- get_map(location = "raleigh")
ggmap(raleigh, extent = "normal")
#The `ggmap` package comes with a "crime" dataset. We will use it here for some plotting
#exercise. Let's only subset violent crimes
violent_crimes <- crime[crime$offense != "auto theft" & crime$offense != "theft" & crime$offense != "burglary", ]
#let's order the violent crimes
violent_crimes$offense <- factor(violent_crimes$offense,
levels = c("robbery", "aggravated assault", "rape", "murder"))
#we can now subset the violent crimes happening within the bounding box seen above
violent_crimes <- violent_crimes[-95.39681 <= violent_crimes$lon & violent_crimes$lon <= -95.34188 &
29.73631 <= violent_crimes$lat & violent_crimes$lat <= 29.78400, ]
#one last step before plotting. We may want to choose a background theme
#among the ones offered by the `ggplot2` package
theme_set(theme_bw(16))
#Finally, using a syntax equal to the `ggplot2` package, let's plot our crimes on top
#of our tile
HoustonMap <- qmap("houston", zoom = 14, color = "bw", legend = "topleft")
HoustonMap +
geom_point(aes(x = lon, y = lat, colour = offense, size = offense),
data = violent_crimes)
## Warning: Removed 12 rows containing missing values (geom_point).
The maps we have created so far are useful to visualize spatial information on top of static basemaps (similar to traditional old-school paper maps). However, because nowadays we are connected to the web all the times, we should consider the added power of designing interactive maps for the web. These typically access the most up-to-date information and can have specialized tools for retrieving information. Moreover, web maps are much more easy to read by the general public that is used to seeing them all the time. A good example of these sort of interactive maps is Google Maps. The main problem most R Users face is that they work on it all the time but do not nenecessarily know how to code in javascript (which would be needed for most web map applications). Luckily, R now offers the opportunity to create nice interactive maps for the web without any knowledge of javascript and HTML! Here we will exlore two main packages: plotGoogleMapsand googleVis.
plotGoogleMapsand googleVis packagesThe plotGoogleMaps package uses the power of Google’s APIs to create intuitive and fully interactive web maps. We can use this package to create stunning HTML pages and later upload them to our websites or share with colleagues. More tutorials and examples similar to these ones can be found HERE and HERE.
We can start by looking at how to plot spatial point data. NOTE: the HTML files generated here, should open automatically in your default browser and saved to a local working directory folder.
#Load point data from a dataset included in the package
data(meuse)
#convert the data frame to a spatial data frame
coordinates(meuse)<-~x+y
#adding Coordinate Referent Sys
proj4string(meuse) <- CRS('+init=epsg:28992')
#Create web map of Point data
m <- plotGoogleMaps(meuse, filename='myMap1.htm')
## [1] "using original PolyCol"
## [1] "using original PolyCol"
### THE HTML FILE WILL OPEN IN YOUR BROWSER AUTOMATICALLY ####
#or you can plot a bubble graph with graduated symbols by zinc values
m <- bubbleGoogleMaps(meuse, zcol='zinc', max.radius = 80, filename='myMap2.htm')
## [1] "using original PolyCol"
### THE HTML FILE WILL OPEN IN YOUR BROWSER AUTOMATICALLY ####
We may be interested in plotting spatial polygon data (choropleth maps). We need to use the R package RColorBrewer to deal with color coding.
#inside the maptools package we can find a shapefile with North Carolina counties
nc <- readShapeSpatial( system.file("shapes/sids.shp",
package="maptools")[1], proj4string=CRS("+proj=longlat
+datum=NAD27"))
#let's make a choropleth map
m <- plotGoogleMaps(nc, zcol="NWBIR74", filename='MyMap3.htm', mapTypeId='TERRAIN',
colPalette= brewer.pal(7,"Reds"),
strokeColor="white")
### THE HTML FILE WILL OPEN IN YOUR BROWSER AUTOMATICALLY ####
Another very useful package to create interactive visualizations, in our case maps, is googleVis. This package lets you can create rich interactive graphics that you can play with locally in Rstudio or in your browser. For more information and examples, please check HERE. The example shown here, is adjusted from the R tutorial for spatial statistics blog. In our example, we are going to import a shapefile with the borders of all the countries in the world found on thematicmapping.org.
library(googleVis)
library(rgdal)
#let's start by reading in the shapefile
polygons <- readOGR("data", layer="TM_WORLD_BORDERS_SIMPL-0.3")
## OGR data source with driver: ESRI Shapefile
## Source: "data", layer: "TM_WORLD_BORDERS_SIMPL-0.3"
## with 246 features
## It has 11 fields
class(polygons)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
#assign this class back to a simple dataframe object
data.poly <- as.data.frame(polygons)
head(data.poly)
## FIPS ISO2 ISO3 UN NAME AREA POP2005 REGION SUBREGION
## 0 AC AG ATG 28 Antigua and Barbuda 44 83039 19 29
## 1 AG DZ DZA 12 Algeria 238174 32854159 2 15
## 2 AJ AZ AZE 31 Azerbaijan 8260 8352021 142 145
## 3 AL AL ALB 8 Albania 2740 3153731 150 39
## 4 AM AM ARM 51 Armenia 2820 3017661 142 145
## 5 AO AO AGO 24 Angola 124670 16095214 2 17
## LON LAT
## 0 -61.783 17.078
## 1 2.632 28.163
## 2 47.395 40.430
## 3 20.068 41.143
## 4 44.563 40.534
## 5 17.544 -12.296
#subset only country name and population fields
data.poly <- data.poly[, c("NAME", "POP2005")]
data.poly$POP2005 <- data.poly$POP2005 / 1000000 #divide population by 1 million
names(data.poly) <- c("Country Name","Population 2005 (millions)")
#now we can use the gvisGeoMap() function to create the desired interactive map
map <- gvisGeoMap(data=data.poly, locationvar = "Country Name",
numvar='Population 2005 (millions)',
options=list(width='800px',heigth='500px',colors="['0x0000ff', '0xff0000']"))
#save map to file using print()
print(map, file="MyMap4.html")
### THE HTML FILE IS SAVED ON YOUR LOCAL FOLDER ####
leafletR packageIn this section, we will show how to use power of the leafletR package to create very nice interactive online maps, without the need of writing a single line of Javascript or HTML5. This package uses the open-source JavaScript library called leaflet. In this sections, you will also start working with GeoJSON files. JSON (JavaScript Object Notation) files are lightweight data-interchange formats, easy for humans to read and write and easy for machines to parse and generateare. They are very commonly used to exchange and save data from web applications. GeoJSON adds encoding features and support for a variety of geographic data structures (i.e. geometries). You can read more about it HERE.
Here we import a polygon boundary of the state of Utah and a set of polygons representing game management units within the state.
deer_unit <- readOGR(dsn = "data", layer = "deerrange_UT") # Game management units
## OGR data source with driver: ESRI Shapefile
## Source: "data", layer: "deerrange_UT"
## with 29 features
## It has 16 fields
#The attribute data is stored in a slot called `@data` so if we just want to look at the data stored in the object:
summary(deer_unit@data)
## FID_ungula UName UNum
## Min. : 1.00 South Slope : 4 Min. : 1.0
## 1st Qu.:15.00 Nine Mile : 2 1st Qu.: 9.0
## Median :24.00 North Slope : 2 Median :14.0
## Mean :24.14 Wasatch Mtns: 2 Mean :15.1
## 3rd Qu.:39.00 West Desert : 2 3rd Qu.:20.0
## Max. :46.00 Book Cliffs : 1 Max. :30.0
## (Other) :16
## SUName SUNum AREA_km2
## none :11 10ab : 1 Min. : 589.1
## Abajo Mountains : 1 11a : 1 1st Qu.: 2655.1
## Anthro : 1 11b : 1 Median : 3968.5
## Bitter Creek, south : 1 12 : 1 Mean : 5873.3
## Bonanza : 1 13a : 1 3rd Qu.: 6637.8
## Currant Creek-Avintaquin: 1 14a : 1 Max. :23325.1
## (Other) :13 (Other):23
## SURF_AREA DEER DEER_SD ELK
## Min. : 180 Min. : 0 Min. : 0.0 Min. : 0
## 1st Qu.: 998 1st Qu.: 1530 1st Qu.: 171.0 1st Qu.: 109
## Median :1781 Median : 8300 Median : 524.0 Median : 1020
## Mean :2004 Mean : 7534 Mean : 704.8 Mean : 1767
## 3rd Qu.:3161 3rd Qu.:11110 3rd Qu.:1089.0 3rd Qu.: 2210
## Max. :4955 Max. :21180 Max. :1875.0 Max. :11180
##
## ELK_SD FAWNS FAWN_SD CALVES
## Min. : 0.0 Min. : 0.00 Min. : 0.000 Min. : 0
## 1st Qu.: 30.0 1st Qu.:54.00 1st Qu.: 5.000 1st Qu.: 0
## Median : 91.0 Median :61.00 Median : 6.300 Median : 0
## Mean :167.4 Mean :54.52 Mean : 6.852 Mean :15
## 3rd Qu.:226.0 3rd Qu.:66.00 3rd Qu.: 9.000 3rd Qu.:42
## Max. :968.0 Max. :79.00 Max. :14.400 Max. :56
##
## CALF_SD ANNRANGE
## Min. :0.000 Min. : 2.562
## 1st Qu.:0.000 1st Qu.: 41.457
## Median :0.000 Median : 178.747
## Mean :1.783 Mean : 551.329
## 3rd Qu.:3.000 3rd Qu.: 639.717
## Max. :8.600 Max. :2742.693
##
names(deer_unit) # The attribute data names
## [1] "FID_ungula" "UName" "UNum" "SUName" "SUNum"
## [6] "AREA_km2" "SURF_AREA" "DEER" "DEER_SD" "ELK"
## [11] "ELK_SD" "FAWNS" "FAWN_SD" "CALVES" "CALF_SD"
## [16] "ANNRANGE"
class(deer_unit)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
#Calculate a new field for vector data, such as the total number of adult deer and elk in each management unit:
deer_unit$tot_deer_elk <- deer_unit$DEER + deer_unit$ELK
summary(deer_unit@data)
## FID_ungula UName UNum
## Min. : 1.00 South Slope : 4 Min. : 1.0
## 1st Qu.:15.00 Nine Mile : 2 1st Qu.: 9.0
## Median :24.00 North Slope : 2 Median :14.0
## Mean :24.14 Wasatch Mtns: 2 Mean :15.1
## 3rd Qu.:39.00 West Desert : 2 3rd Qu.:20.0
## Max. :46.00 Book Cliffs : 1 Max. :30.0
## (Other) :16
## SUName SUNum AREA_km2
## none :11 10ab : 1 Min. : 589.1
## Abajo Mountains : 1 11a : 1 1st Qu.: 2655.1
## Anthro : 1 11b : 1 Median : 3968.5
## Bitter Creek, south : 1 12 : 1 Mean : 5873.3
## Bonanza : 1 13a : 1 3rd Qu.: 6637.8
## Currant Creek-Avintaquin: 1 14a : 1 Max. :23325.1
## (Other) :13 (Other):23
## SURF_AREA DEER DEER_SD ELK
## Min. : 180 Min. : 0 Min. : 0.0 Min. : 0
## 1st Qu.: 998 1st Qu.: 1530 1st Qu.: 171.0 1st Qu.: 109
## Median :1781 Median : 8300 Median : 524.0 Median : 1020
## Mean :2004 Mean : 7534 Mean : 704.8 Mean : 1767
## 3rd Qu.:3161 3rd Qu.:11110 3rd Qu.:1089.0 3rd Qu.: 2210
## Max. :4955 Max. :21180 Max. :1875.0 Max. :11180
##
## ELK_SD FAWNS FAWN_SD CALVES
## Min. : 0.0 Min. : 0.00 Min. : 0.000 Min. : 0
## 1st Qu.: 30.0 1st Qu.:54.00 1st Qu.: 5.000 1st Qu.: 0
## Median : 91.0 Median :61.00 Median : 6.300 Median : 0
## Mean :167.4 Mean :54.52 Mean : 6.852 Mean :15
## 3rd Qu.:226.0 3rd Qu.:66.00 3rd Qu.: 9.000 3rd Qu.:42
## Max. :968.0 Max. :79.00 Max. :14.400 Max. :56
##
## CALF_SD ANNRANGE tot_deer_elk
## Min. :0.000 Min. : 2.562 Min. : 0
## 1st Qu.:0.000 1st Qu.: 41.457 1st Qu.: 2604
## Median :0.000 Median : 178.747 Median : 9786
## Mean :1.783 Mean : 551.329 Mean : 9302
## 3rd Qu.:3.000 3rd Qu.: 639.717 3rd Qu.:13503
## Max. :8.600 Max. :2742.693 Max. :32360
##
Many shapefiles have significant detail that results in very large GeoJSON (JavaScript Object Notation) files. If you do not simplify the topology, a file can end up being several MB. Simplifying a topology represents an important GIS technique to reduce the vertex count in features that were captured in too much detail. More detail can be found HERE. Here, we will use the rgeos package to simplify topology. NOTE: if preserving topology is crucial, consider using tools in QGIS or GRASS GIS. For this demonstration we will restrict ourselves to R
library(rgeos)
#save the data slot
deer_unit_sub <- deer_unit@data[, c("UName", "UNum", "AREA_km2", "tot_deer_elk")]
#simplification yields a SpatialPolygons class
deer_unit <- gSimplify(deer_unit, tol=0.01, topologyPreserve=TRUE)
class(deer_unit)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
#to write to geojson we need a SpatialPolygonsDataFrame
deer_unit <- SpatialPolygonsDataFrame(deer_unit, data=deer_unit_sub)
class(deer_unit)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
head(deer_unit@data)
## UName UNum AREA_km2 tot_deer_elk
## 0 Zion 29 4454.250 9786
## 1 Paunsaugunt 27 3873.350 5794
## 2 Kaiparowitz 26 8127.691 550
## 3 Pine Valley 30 6753.223 13210
## 4 Panguitch Lake 28 2286.769 10035
## 5 Mt. Dutton 24 1702.916 3870
We are now ready to create the GeoJSON file. We will also create the variable intervals (cuts) as we want in order to map it.
#write data to GeoJSON
writeOGR(deer_unit, dsn="./data/DeerElkGeoJson", layer="DeerElkGeoJson", driver="GeoJSON")
#a GeoJSON datasource is translated to single OGRLayer object with pre-defined name OGRGeoJSON"
ogrInfo("./data/DeerElkGeoJson", "OGRGeoJSON")
## Source: "./data/DeerElkGeoJson", layer: "OGRGeoJSON"
## Driver: GeoJSON; number of rows: 29
## Feature type: wkbPolygon with 2 dimensions
## Extent: (-114.0572 37.00046) - (-109.0348 42.00169)
## CRS: +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs
## Number of fields: 4
## name type length typeName
## 1 UName 4 0 String
## 2 UNum 0 0 Integer
## 3 AREA_km2 2 0 Real
## 4 tot_deer_elk 2 0 Real
#Divide the variable of interest into bins using custom cuts (in this case quintiles)
cuts <- round(quantile(deer_unit$tot_deer_elk, probs = seq(0, 1, 0.20), na.rm = FALSE), 0)
We are almost ready to create the final interactive web map. First, we want to choose a number of fields that we want displayed in a pop-up balloon inside our map. We may also want to define a graduated color style to plot our attribute of interest. For that, we will use the styleGrad() function of the leafletR package.
#Choose fields to include in the popup
popup <- c("UName", "UNum", "AREA_km2", "tot_deer_elk")
#Graduated style symbology based on an attribute
sty <- styleGrad(prop="tot_deer_elk", breaks=cuts, closure="left", style.par="col",
style.val=rev(heat.colors(5)), leg="Deer & Elk Population", lwd=1)
#Create the map and load into browser
map <- leaflet(data="./data/DeerElkGeoJson", style=sty, title="DeerElk",
base.map="osm", incl.data=TRUE, popup=popup)
#To look at the map you can either find the .html file in your local folder, or call it using
#browseURL(map)
In this module, we will work on a dataset containing abundance records of California bay laurel (UMCA) trees in the Sonoma county, California. After reading in and cleaning the dataset, we will create a leaflet heat map using a JavaScript plugin. You can find plenty of more examples and documentation HERE.
#read in data table
tree_dat <- read.csv("./data/plots_umca_infection.csv", header=T)
#read plot locations shapefile
plots <- readOGR(dsn='./data', layer='plot_202_latlon')
## OGR data source with driver: ESRI Shapefile
## Source: "./data", layer: "plot_202_latlon"
## with 202 features
## It has 5 fields
#let's subset data for year 2012 only and 3 variables of interest.
tree_dat_2012 <- tree_dat[tree_dat$year == 2012, c("plot", "year", "tot_bay")]
head(tree_dat_2012)
## plot year tot_bay
## 9 ANN01 2012 26
## 19 ANN02 2012 10
## 29 ANN03 2012 45
## 39 ANN04 2012 13
## 49 ANN05 2012 14
## 59 ANN06 2012 10
We now need to add the geographic coordinates of each plot. We have 202 plot locations in the shapefile but our dataset only has 179. We can use the R function merge() to keep only the matching ones.
tree_dat_2012 <- merge(tree_dat_2012, plots@data, by.x="plot", by.y="PLOT_ID")
head(tree_dat_2012)
## plot year tot_bay X Y POINT_X POINT_Y
## 1 ANN01 2012 26 532568 4255157 -122.6268 38.44411
## 2 ANN02 2012 10 532557 4254989 -122.6269 38.44259
## 3 ANN03 2012 45 533352 4254736 -122.6178 38.44029
## 4 ANN04 2012 13 532314 4254707 -122.6297 38.44006
## 5 ANN05 2012 14 533109 4254426 -122.6206 38.43750
## 6 ANN06 2012 10 531353 4254410 -122.6408 38.43742
#remove unnecessary variables
tree_dat_2012 <- tree_dat_2012[ , c("POINT_Y","POINT_X","tot_bay")]
#rename the coordinates fields to lat and lon
names(tree_dat_2012)[1:2] <- c("lat", "lon")
head(tree_dat_2012)
## lat lon tot_bay
## 1 38.44411 -122.6268 26
## 2 38.44259 -122.6269 10
## 3 38.44029 -122.6178 45
## 4 38.44006 -122.6297 13
## 5 38.43750 -122.6206 14
## 6 38.43742 -122.6408 10
Now, we can create a new Leaflet map instance and set the plot extent based on our coordinates centroid.
#create a new leaflet map instance
Lmap <- Leaflet$new()
#set the view and zoom to the desired study area. Let's center it on our mean lat-lon coordinates
Lmap$setView(c(mean(tree_dat_2012$lat), mean(tree_dat_2012$lon)), zoom=10)
#add a terrain basemap
Lmap$tileLayer(provider = "MapQuestOpen.OSM")
Before generating the interactive web map, we need to convert our dataset to JSON format using the toJSONArray2() function.
tree_dat <- toJSONArray2(na.omit(tree_dat_2012), json = F, names = F)
#let's print out the first two elements of the JSON file
cat(toJSON(tree_dat[1:2]), '\n')
## [[38.444108098,-122.626806828,26],[38.4425944789,-122.626940672,10]]
The heat map is created thank to a leaflet-heat plugin, wrote by Vladimir Agafonkin. We only need to point to the plugin web address to use it with our map.
#add leaflet-heat plugin. Thanks to Vladimir Agafonkin
Lmap$addAssets(jshead = c(
"http://leaflet.github.io/Leaflet.heat/dist/leaflet-heat.js"
))
#add javascript to modify underlying chart
Lmap$setTemplate(afterScript = sprintf("
<script>
var addressPoints = %s
var heat = L.heatLayer(addressPoints).addTo(map)
</script>
", rjson::toJSON(tree_dat)
))
Lmap
## <iframe src=' ./data/figures/leaflet9-1.html ' scrolling='no' frameBorder='0' seamless class='rChart leaflet ' id=iframe- chartd8c7649c448 ></iframe> <style>iframe.rChart{ width: 100%; height: 400px;}</style>
### CHECK YOUR BROWSER INSIDE THE figures FOLDER TO SEE THE MAP ####
The following tutorial is based on the rMaps online tutorial used to create animated choropleths. The data selected are different from the original. The first step to creating any visualization is getting the data. Let us fetch time-series data on burglaries in the US, from [Quandl] (https://www.quandl.com/).
#let's download the table referring to burglaries
rbData = Quandl("FBI_UCR/USCRIME_TYPE_BURGLARIES")
rbData[1:10, 1:5]
## Year Alabama Alaska Arizona Arkansas
## 1 2010-12-31 42034 3105 50771 32511
## 2 2009-12-31 48844 3600 54308 34753
## 3 2008-12-31 50379 3232 58125 33772
## 4 2007-12-31 45360 3683 59290 31874
## 5 2006-12-31 44561 4122 58768 32450
## 6 2005-12-31 43473 4131 56328 30118
## 7 2004-12-31 44666 3773 56885 30151
## 8 2003-12-31 43245 3874 58613 25016
## 9 2002-12-31 42578 3908 59087 23229
## 10 2001-12-31 40642 3847 54821 22196
The dataset is in the wide-form. The first step is to convert it into the long-form (more convenient for visualization purposes). Moreover, we remove data for the US as a whole, as well as for DC, so that the crime rates across entities (states) are comparable
#let's use the melt() function to restructure the data in wide form
datm <- melt(rbData, 'Year',
variable.name = 'State',
value.name = 'Crime'
)
datm <- subset(na.omit(datm),
!(State %in% c("United States", "District of Columbia"))
)
head(datm)
## Year State Crime
## 1 2010-12-31 Alabama 42034
## 2 2009-12-31 Alabama 48844
## 3 2008-12-31 Alabama 50379
## 4 2007-12-31 Alabama 45360
## 5 2006-12-31 Alabama 44561
## 6 2005-12-31 Alabama 43473
Crime rates need to be categorized (discrete classes). One way to do this is to divide them into quartiles.
#let's use the function transform() to change our data frame
datm2 <- transform(datm,
State = state.abb[match(as.character(State), state.name)],
fillKey = cut(Crime, quantile(Crime, seq(0, 1, 1/4)), labels = LETTERS[1:4]),
Year = as.numeric(substr(Year, 1, 4))
)
Each quartile needs to be associated with a fill color chosen from a palette. We can use the RColorBrewer package to deal with custom color fill.
fills = setNames(
c(brewer.pal(4, 'OrRd'), 'white'),
c(LETTERS[1:4], 'defaultFill')
)
We need to convert the data frame into a list of lists. We will use Hadley’s plyr package to simplify the code.
dat2 <- dlply(na.omit(datm2), "Year", function(x){
y = toJSONArray(x, json = F)
names(y) = lapply(y, '[[', 'State')
return(y)
})
names(dat2)
## [1] "1960" "1961" "1962" "1963" "1964" "1965" "1966" "1967" "1968" "1969"
## [11] "1970" "1971" "1972" "1973" "1974" "1975" "1976" "1977" "1978" "1979"
## [21] "1980" "1981" "1982" "1983" "1984" "1985" "1986" "1987" "1988" "1989"
## [31] "1990" "1991" "1992" "1993" "1994" "1995" "1996" "1997" "1998" "1999"
## [41] "2000" "2001" "2002" "2003" "2004" "2005" "2006" "2007" "2008" "2009"
## [51] "2010"
#dat2[["1960"]] to inspect a list element
We can now create a simple choropleth map of crime rates for a given year. The Datamaps reference class gives us simple bindings to the DataMaps library
options(rcharts.cdn = TRUE)
map <- Datamaps$new()
map$set(
dom = 'chart_1',
scope = 'usa',
fills = fills,
data = dat2[["1980"]],
legend = TRUE,
labels = TRUE
)
map
## <iframe src=' ./data/figures/choro6-1.html ' scrolling='no' frameBorder='0' seamless class='rChart datamaps ' id=iframe- chart_1 ></iframe> <style>iframe.rChart{ width: 100%; height: 400px;}</style>
Use a customized wrapper function that absorbs js code to produce a dynamic choropeth map.
map2 <- ichoropleth(Crime ~ State,
data = datm2[,1:3],
pal = 'OrRd', #color ramp
ncuts = 4, #quartiles
animate = 'Year'
)
map2
## <iframe src=' ./data/figures/choro7-1.html ' scrolling='no' frameBorder='0' seamless class='rChart datamaps ' id=iframe- chartd8c5be46995 ></iframe> <style>iframe.rChart{ width: 100%; height: 400px;}</style>